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Abstract. Moduli of rings and quadrilaterals are frequently applied in geomet- 
ric function theory, see e.g. the Handbook by Kiihnau. Yet their exact values 
are known only in a few special cases. Previously, the class of planar domains 
with polygonal boundary has been studied by many authors from the point of 
view of numerical computation. We present here a new /ip-FEM algorithm for the 
computation of moduli of rings and quadrilaterals and compare its accuracy and 
performance with previously known methods such as the Schwarz-Christoffel Tool- 
box of Driscoll and Trefethen. We also demonstrate that the /ip-FEM algorithm 
applies to the case of non-polygonal boundary and report results with concrete 
error bounds. 
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1. Introduction 

Plane domains with piecewise-smooth boundary curves occur in applications to 
electronics circuit design, airfoil modelling in computational fluid dynamics, com- 
puter vision and various other problems of engineering and science [23l [271 EH [Ml 
[371 [32] • We assume that the domain is bounded and that there are either one or two 
simple (and nonintersecting) boundary curves. The domain is then either simply or 
doubly connected. For the mathematical modelling of these domains it is usually 
convenient to map the domains conformally onto "canonical domains" as simple as 
possible: the unit disk D = e C : 1^1 < 1} or the annulus {z E C : r < \z\ < 1} . 
Sometimes a rectangle is preferable to the unit disk as a canonical domain. The 
existence of these canonical conformal mappings is guaranteed by classical results of 
geometric function theory but the construction of this mapping in a concrete appli- 
cation case is usually impossible. Therefore one has to resort to numerical conformal 
mapping methods for which there exists an extensive literature flEi 128] [311 137] . The 
Schwarz-Christoffel (SC) Toolbox of Driscoll [l7j, based on the software of Trefethen 
|41j . is in wide use for numerical conformal mapping applications. 

In the doubly connected case, one might be interested in only knowing the inner 
radius r of the canonical annulus. For instance this occurs if we wish to compute the 
electric resistance of a ring condenser. In this situation the conformal mapping itself 
is not needed if we are able to flnd the inner radius r by some other method. It is a 
classical fact that the inner radius r can be obtained in terms of the solution of the 
Dirichlet problem for the Laplace equation in the original domain with the boundary 
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value on one boundary component and the boundary value 1 on the other one. 
This fact is just one way of saying that the modulus of a ring domain is conformally 
invariant: for the canonical annulus {z ^ C : r < \z\ < 1} the modulus is equal to 
log(l/r) . This idea reduces the problem of computing the number r to the problem 
of numerical approximation of the solutions of Laplace equation in ring domains. In 
the paper [9] this method was applied to several concrete examples of ring domains 
for which numerical results were reported, too. Again, it is also possible to use the 
Schwarz-Christoffel method for doubly connected domains [23]. 

We next consider the case of simply connected plane domains. For such a domain 
D and for a quadruple {zi, Z2, z^, z^} of its boundary points we call {D; zi, Z2, z^, z^) a 
quadrilateral if zi, 2:2, Zs, z^ occur in this order when the boundary curve is traversed 
in the positive direction. The points z^^k = 1, ..,4, are called the vertices and the 
part of the oriented boundary between two successive vertices such as Zi and Z2 is 
called a boundary arc (21, Z2) ■ The modulus M{D; zi, Z2, Z3, Z4) of the quadrilateral 
{D; zi, Z2, Z3, Z4) is defined to be the unique /i > for which there exists a conformal 
mapping of D onto the rectangle with vertices 1 + ih, ih, 0, 1 such that the points 
Zi, Z2, Z3, Z4 correspond to the vertices in this order. This conformal mapping is called 
the canonical conformal mapping associated with the quadrilateral. As in the case of 
doubly connected domains discussed above, it is well-known that the computation of 
the modulus h of the quadrilateral may be reduced to solving the Dirichlet-Neumann 
boundary value problem in the original domain D with the Dirichlet boundary values 
1 on the boundary arc {zi, Z2) and on the arc {z^, Z4) and Neumann boundary values 
on the arcs (2:3, z^) and (2:4, zi) . 

Conformal moduli of rings and quadrilaterals have independent theoretical in- 
terest because of their crucial role in the theory of quasiconformal mappings |29] . 
These quantities are closely related to certain physical constants, e.g. they play an 
important role in applications involving the measurement of resistance values of inte- 
grated circuit networks. But the problem of computing the moduli is also interesting 
in the wider engineering context. The reciprocal identities (4.1 ) and (6.3 ) can be used 
to generate test cases for curvilinear Lipschitz domains and thus should be standard 
tools in the FEM-software development community. Unfortunately these identities 
are missing from the introductory FEM textbooks. Even though our interest lies in 
the high-order methods, these test cases are equally valid for any numerical PDE 
methods and mesh adaptation in particular. 

One specific application area of the algorithms presented here is the simulation 
of measurements for the 2D electrical impedance tomography (EIT) [25]. In EIT 
problems a number of electrodes are placed on the boundary of the domain and 
current patterns are considered between every pair of them. Indeed, computing the 
moduli can be considered as a very crude model for the so-called EIT background 
forward problem. In general, the meshes for the EIT forward problems can be adapted 
using the approaches outlined below. High level of accuracy is necessary for precise 
control of artificial noise in the simulations. 

A general observation about the literature seems to be that reported numerical 
values of the moduli of concrete quadrilaterals (or ring domains) are hard to find 
in the literature. Perhaps the longest list of numerical results is given in [9j where 
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pointers to earlier literature may be found. The recent book [M] lists also many such 
numerical values. In our opinion a catalogue of these numerical values in the simplest 
cases would be desirable for instance for reference purposes. The book [37] and the 
paper [3^ p. 127] list certain engineering formulas which have been applied in VLSI 
circuit design. 

An outline of the structure of this paper now follows. First, in Section 2 we 
describe the methods used in this paper. In Section 3 we discuss in detail the various 
FEM methods used here, in particular the hp-method which was implemented and 
applied to generate some of the results reported below. Another method we use is the 
/i-adaptive software package AFEM of K. Samuelsson, which implements an adaptive 
FEM method and which was previously used in pj. In the present paper we use the 
AFEM method to compute the modulus of a quadrilateral whereas in |UJ it was 
used merely for the computation of the moduli of ring domains. In Section 4 a test 
problem for quadrilaterals is described together with its analytic solution, following 
|22j . This analytic solution requires, however, an application of a numerical root 
finding program. Accordingly, this formula is analytic-numeric in its character. In 
Section 5 we check several methods against this analytic formula in a test involving a 
family of convex quadrilaterals. The methods discussed are the analytic formula from 
|22j . the Schwarz-Christoffel Toolbox of pTlfTS] , the AFEM method of Samuelsson [9] 
and the present hp-method. On the basis of these experiments, an accuracy ranking 
of the methods is given in Section 5. In Section 6 the more general case of polygonal 
quadrilaterals is investigated, in particular L-shaped domains, and the results are 
compared to the literature. In Section 7 we discuss the computation of the modulus 
of a ring domain in a few special cases. For instance, for "the cross in square" 
ring domain considered previously in [9l Example 4] we now obtain much improved 
accuracy. In Section 8 we compute some examples with the hp-FEM which are 
difficult for other methods. In Section 9 our results and discoveries are summarized. 

2. Methods 

The following problem is known as the Dirichlet- Neumann problem. Let D be a 
region in the complex plane whose boundary dD consists of a finite number of regular 
Jordan curves, so that at every point, except possibly at finitely many points, of the 
boundary a normal is defined. Let dD = AUB where A, B both are unions of Jordan 
arcs. Let ipA, i^B be a real-valued continuous functions defined on A, B, respectively. 
Find a function u satisfying the following conditions: 

(1) u is continuous and differentiable in D. 

(2) u{t) = ipAit), for all t e A. 

(3) If d/dn denotes differentiation in the direction of the exterior normal, then 

d 

—u(t) = tpBit), for all t e B. 
on 

2.1. Modulus of a quadrilateral and Dirichlet integrals. One can express the 
modulus of a quadrilateral {D; zi, Z2, z^, z^) in terms of the solution of the Dirichlet- 
Neumann problem as follows. Let ■yj,j = 1, 2, 3, 4 be the arcs of dD between (^4, zi) , 
{zi,Z2) , (-22,-23) , (-23,^4), respectively. If u is the (unique) harmonic solution of the 
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Dirichlet-Neumann problem with boundary values of u equal to on 72, equal to 1 
on 74 and with du/dn = on 71 U 73 , then by ll, p. 65/Thm 4.5]: 



2.3. Modulus of a ring domain and Dirichlet integrals. Let E and F be two 

disjoint compact sets in the extended complex plane Cqo- Then one of the sets E, 
F is bounded and without loss of generality we may assume that it is . If both E 
and F are connected and the set R = Coo \{E U F) is connected, then R is called a 
ring domain. In this case i? is a doubly connected plane domain. The capacity of R 
is defined by 



where the infimum is taken over all nonnegative, piecewise differentiable functions u 
with compact support in RU E such that u = 1 on E. It is well-known that the har- 
monic function on R with boundary values 1 on E and on F is the unique function 
that minimizes the above integral. In other words, the minimizer may be found by 
solving the Dirichlet problem for the Laplace equation in R with boundary values 1 
on the bounded boundary component E and on the other boundary component F . 
A ring domain R can be mapped conformally onto the annulus {z : < \z\ < 1}, 
where M = M{R) is the conformal modulus of the ring domain R. The modulus and 
capacity of a ring domain are connected by the simple identity M{R) = 27r/cap R. For 
more information on the modulus of a ring domain and its applications in complex 
analysis the reader is referred to [H [231 EH EHJ El] . 

In |3H Chapter 3] N. Papamichael and N. Stylianopoulos describe the so-called 
domain decomposition method for the computation of the modulus of a quadrilat- 
eral which is designed for the case of elongated quadrilaterals and applies e.g. to 
polygonal quadrilaterals that can be decomposed into simple pieces whose moduli 
can be estimated. As an example they consider a spiraling quadrilateral that can be 
decomposed into a "sum" of 13 trapezoids and report results that are expected to 
be correct up to 7 decimal places. Therefore, this method seems very attractive for 
the computation of the modulus of a special class of quadrilaterals. A key feature 
of the method is that it reduces the numerical difficulties caused by the crowding 
phenomenon for this special class of quadrilaterals. 

2.4. Classification of methods for numerical computing. For the computation 
of the modulus of a quadrilateral or of a ring domain there are two natural approaches 

(1) methods based on the definition of the modulus and on the use of a conformal 
mapping onto a canonical rectangle or annulus, 

(2) methods that give only the modulus, not the canonical conformal map. 

In some sense, methods of class (1) give a lot of extra information, namely the 
conformal mapping - all we want is a single real number. Methods of class (2) rely 
on solving the Dirichlet-Neumann boundary value problem or Dirichlet problem for 
the Laplace equation as described above. 

In this paper we will mainly use methods of type (2) that make use of adaptive 
FEM methods for solving the Laplace equation. 



(2.2) 
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2.5. Review of the literature on numerical conformal mapping. With the ex- 
ception of a few special cases, both of the above methods lead to extensive numerical 
computation. For both classes of methods there are several options in the literature, 
see for instance the bibliography of [0] . Various aspects of the theory and practice of 
numerical conformal mapping are reviewed in the monographs [IHl EHl EH EZ] • See 
also the authoritative surveys [201 [82lli2l HH] . 

Recently numerical conformal mappings have been studied from various points 
of view and in various applications by many authors, see e.g. [21 El |T3l HH |15l [26l 



In the paper [9] the modulus of a ring domain was computed with the help of 
the software package AFEM of K. Samuelsson, based on an /i-adaptive finite element 
method. It can be easily applied to compute the modulus of a quadrilateral. 

In this section we describe the high-order p-, and hp-finite element methods and 
report the results of numerical computation of the moduli of a number of quadrilat- 
erals. The paper of Babuska and Suri [7| gives an accessible overview of the method, 
for a more detailed exposition we refer to Schwab [38], and for those familiar with 
engineering approach the book by Szabo and Babuska is recommended. 

In the /i-version or standard finite element method, the unknowns or degrees 
of freedom are associated with values at specified locations of the discretization of 
the computational domain, that is, the nodes of the mesh. In the p-method, the 
unknowns are coefficients of some polynomials that are associated with topological 
entities of the elements, nodes, sides, and interior. Thus, in addition to increasing 
accuracy through refining the mesh, we have an additional refinement parameter, the 
polynomial degree p. 

Let us next define a p-type quadrilateral element. The construction of triangles 
is similar and can be found from the references given above. 

3.1. Shape functions. Many different selections of shape functions are possible. We 
follow Szabo and Babuska [IQ] and present the so-called hierarchic shape functions. 

Legendre polynomials of degree n can be defined using a recursion formula 



For our purposes the central polynomials are the integrated Legendre polyno- 
mials 



EIlESlEe]. 



3. p-, AND /ip- FINITE ELEMENT METHOD 



(3.2) (n + l)P„,+i(x) - (2n + l)xP„(x) + nP„_i(x) = 0, Po(x) = 1. 
The derivatives can similarly be computed using a recursion 

(3.3) (1 - a;2)p;(a;) = -nxP^ix) + nP„_i(x). 



(3.4) 
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The normalizing coefficients are cfiosen so tfiat 

(3.6) L^^^^^'^- 

We can now define tfie sfiape functions for a quadrilateral reference element. 
The shape functions are divided into three categories: nodal shape functions, side 
modes, and internal modes. 

3.7. Nodal shape functions. There are four nodal shape functions. 



iVi(e,^) = 






^(1 + 0(1-^), 




^(1 + 0(1 + ^7), 




|(1-0(1 + ^)- 



Taken alone, these shapes define the standard four-node quadrilateral finite element. 

3.8. Side shape functions. There are 4(p — 1) side modes associated with the sides 
of a quadrilateral (p > 2). 





1 

2' 




t = 2,.. 


■ ■ ,P, 




1 

2' 


[i + OUv), 


i = 2,.. 


■■,P, 




1 

2' 




i = 2,. 


■■,P, 




1 

2' 




t = 2,.. 


■ ,P- 



3.9. Internal shape functions. For the internal modes we have two options. The 

so-called trunk space has (p — 2)(p — 3)/2 shapes 

(3.10) N',{^,r]) = UOMv), hJ>2, t + j=4,5,...,p, 
whereas the full space has {p — l)(p — 1) shapes 

(3.11) N^^{^,r]) = UOMv), ^ = 2,...,p, j = 2,...,p, 

where in both cases the index k depends on the chosen convention. In this paper 
we always use the full space. The internal shape functions are often referred to as 
bubble-functions . 

3.12. Parity problem. The Legendre polynomials have the property Pn{—x) = 
{—l)^Pn{x). In 2D all internal edges of the mesh are shared by two different ele- 
ments. We must ensure that each edge has the same global parameterization in both 
elements. This additional book-keeping is not necessary in the standard h-FEM. 
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3.13. Resource requirements. We have seen that the number of unknowns in a 
p-tjpe quadrilateral is {p + 2){p + 3)/2 or Ap + {p — 1)^ if the internal modes are 
from trunk or full space, respectively. To compensate this, the number of elements 
is naturally taken to be as small as possible. Indeed, when the mesh is adapted in 
a suitable way, the dimension of the overall linear system can be significantly lower 
than in the corresponding /i-method. However, the matrices tend to be denser in the 
p-method, so the space requirements in relation to the dimension of the linear system 
are greater for the p-method. 

3.14. Proper grading of the meshes. For a certain class of problems it can be 
shown that if the mesh and the elemental degrees have been set optimally, we can 
obtain exponential convergence. A geometric mesh is such that each successive layer 
of elements changes in size with some geometric scaling factor a, toward some point 
of interest. In this case, the points of interest are always corner points. 

The following theorem is due to Babuska and Guo [5]. Note that construction 
of appropriate spaces is technical. For rigorous treatment of the theory involved see 
Schwab [38], Babuska and Guo [Ij and references therein. 

3.15. Theorem. Let Q G M."^ be a polygon, v the FEM-solution, and let the weak 
solution Uq be in a suitable countably normed space where the derivatives of arbitrarily 
high order are controlled. Then 



where C and b are independent of N, the number of degrees of freedom. Here v is 
computed on a proper geometric mesh, where the orders of individual elements depend 
on their originating layer, such that the highest layers have the smallest orders. 
The result also holds for constant polynomial degree distribution. 

Let us denote the number of the highest layer with u, the nesting level. Using 
this notation we can refer to geometric meshes as {a, i/)-meshes. 

In Figure [T] we show a geometric mesh template for a non-convex quadrilateral. 
Here we require that each node lies at the end point of an edge and that the meshlines 
follow the guidelines of the geometric meshes. 

In Figure [2] a sequence of real p-type meshes is shown. The template mesh serves 
also as a pure p-type mesh where the approximation properties are changed only by 
varying the polynomial degree. In the middle and rightmost meshes the number of 
elements is the same because the nesting level is the same, only the scaling factor 
changes. 

3.16. Generating geometric meshes. Here we consider generation of geometric 
meshes in polygonal domains. We use the following two-phase algorithm: 

(1) Generate a minimal mesh (triangulation) where the corners are isolated with 

a fixed number of triangles depending on the interior angle, 6 so that the 

refinements can be carried out independently: 

• 6 < 71/2: one triangle, 

• IT /2 < 6 < 7i: two triangles, and 

• 71 < 6: three triangles. 
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Figure 1 . Geometric mesh for a general quadrilateral. 




Figure 2. Graded meshes: Effect of the scaling factor. From left to 
right, template mesh, (a, z/) = (1/2,3), (a, i/) = (1/6,3). 




Figure 3. Three sample meshes used in numerical experiments below. 
Note the refinement of the mesh structure close to the corner points. 

(2) Every triangle attached to a corner is replaced by a refinement, where the 
edges incident to the corner are split as specified by the scaling factor a. This 
process is repeated recursively until the desired nesting level v is reached. 
Note that the mesh may include quadrilaterals after refinement. 

In Figure [2] we can also see our preferred element subdivisions: triangle to 
(quadrilateral, triangle)-pair, and quadrilateral to three quadrilaterals. These two 
rules are sufficients for our purposes since we always grade toward a corner point. 
Using this, we can derive a simple estimate for the number of degrees of freedom iV. 
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Figure 4. Final geometric or (0.15, 12)-mesh. Due to small a only 
first two levels are visible. 




Figure 5. Curved boundary mapping. 



Letting T denote the number of elements in the initial mesh, and C the number of 
corners in the domain (or those used in refining): 

(3.17) iV ~ (T + 6Cz/)p2^ 

where the constant 6 is the product of the maximal number of elements surrounding 
a corner, 3, and the maximal number of new elements per level, 2. 

Finally, in Figure [3] three minimal meshes and in Figure |4] one final mesh are 
shown. 

3.18. Domains with curved boundaries. Since we want to use as large elements 
as possible, it is important to represent curved boundary segments accurately. The 
linear blending function method of Gordon and Hall [2T] is our choice for this purpose. 

In the general case all sides of an element can be curved as in Figure [5] We 
assume that every side is parameterized: 

(3.19) X = Xi(t), y = yi{t), —l<t<l, i = 1, . . . , number of sides 
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Using capital letters as coordinates of the corner points, (Xj,Fj), we can write the 
mapping for the global x-coordinates of a quadrilateral as 

a; = ^(1 - ^)^i(0 + ^(1 + + + + ^(1 - O^M 

(3.20) _ 1(1 _ ^)(i _ ^)Xi - 1(1 + 0(1 + V)X2 - J(l + 0(1 + V)X3 

-^(l-0(l + r])X4, 

and symmetrically for the y-coordinate. Note, that if the side parameterizations 
represent straight edges, the mapping simplifies to the standard bilinear mapping of 
quadrilaterals. 

In the following we always use exact representation of the geometry which im- 
plies that in the ensuing mesh grading process no approximation of geometry is 
necessary. Here the mesh generation of the curved domains is template-based, thus 
the changes in curvature are not automatically dealt with. For a highly accessible 
review of the p-method mesh generation issues we refer to [30] . 



4. Convex quadrilateral 

In this section our goal is to introduce a test problem, whose solution is de- 
termined by a transcendental equation. This equation can be numerically solved 
to the desired accuracy and we will use this to check the validity of the numerical 
methods we use as well as to obtain an experimental estimate for their accuracy. 
The test problems we consider are convex polygonal quadrilaterals. The simplest 
such quadrilateral consists of the four vertices and the line segments joining the ver- 
tices. Let zi, Z2, Z3, Z4 G C be distinct points and suppose that the polygonal line 
that results from connecting these points by segments in the order 2:1,^2,^3,-^4,-^1 
forms the positively oriented boundary of a domain Q. For simplicity, we denote by 
QM{zi, Z2, Zs, Z/i) the modulus M{Q; zi, Z2, z^, z^). Then the modulus is a conformal 
invariant in the following sense: If / : Q — ?• fQ is a conformal mapping onto a Jordan 
domain fQ, then / has a homeomorphic extension to the closure Q (also denoted by 
/) and 

M(g; 21, 2:2,^3, ^4) = M{fQ-J{zi)J{z2),fiz3)J{zi)). 
It is clear by geometry that the following reciprocal identity holds: 

(4.1) M{Q] zi, Z2, Z3, Z4^)M{Q; Z2, Z3, Z4, zi) = 1. 

If h: C — )■ C is a translation, rotation, or stretching, then the piecewise linear nature 
of the boundary is preserved and we can write the conformal invariance simply as 

QM(^i, Z2, ^3, -24) = QM(/(;2i), f{z2), f{z,), f{z,)) . 

There are two particular cases, where we can immediately give QM(zi, Z2, -23, -24). 
The first cases occurs, when all the sides are of equal length (i.e. the quadrilat- 
eral is a rhombus) and in this case the modulus is 1 , see [22] . In the second case 
{Q; zi, Z2, -23, -24) is the rectangle {Q; 1+ih, ih, 0, 1), /i > 0, and QM{l+ih, ih, 0, 1) = h. 
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4.2. Basic identity. In 2.11] some identities satisfied by the function QM(a, b, 0, 1) 
were pointed out. We will need here the following one, which is the basic reciprocal 
identity (4.1) rewritten for the expression QM : 

(4.3) QM(a, b, 0, 1) ■ QM((6 - l)/(a - 1), 1/(1 - a), 0, 1) = 1 . 

We shall consider here the following particular cases of this reciprocal identity: 
(a) parallelogram, (b) trapezoid with angles (7r/4, 37r/4, 71/2, 7r/2), and (c) a convex 
polygonal quadrilateral. Note that for the cases (a) and (b) the formula is less 
complex than for the general case (c). 

4.4. The hypergeometric function and complete elliptic integrals. Given 
complex numbers a,b, and c with c ^ 0,-1,-2,..., the Gaussian hypergeometric 
function is the analytic continuation to the slit plane C \ [1, oo) of the series 



4.5 Fa,&;c;z =2^1 a,&;c;z =E f ^ 

^-^ (c, n) n! 

n=0 ^ ' ^ 

Here (a, 0) = 1 for a ^ 0, and (a, n) is the shifted factorial function or the Appell 
symbol 

(a, n) = a{a + l)(a + 2) ■ ■ ■ (a + n — 1) 

for n E N \ {0}, where N = {0,1,2,...} and the elliptic integrals %{r),%'{r) are 
defined by 

X(r) = 1^(1/2, 1/2; 1; r^), X'(r) = X(r'), and r' = Vl-r^. 

Some basic properties of these functions can be found in [1]. 
4.6. Parallelogram. For t E (0, vr) and h > let 

g{t,h) = QM{l + he'\he'\0,l). 
An analytic expression for this function has been given in [3l 2.3]: 

(4.7) g{t,h) = X'{rt/^)/X{rt/^), 
where 

(4.8) r„ = /i-i \ , for < a < 1/2, 



2 sin(7ra) 

and the decreasing homeomorphism jia'- (0, 1) — )■ (0, oo) is defined by 

TT F(a, 1 — a; 1; 1 — r^) 



(4.9) /i«(r) 



2sin(7ra) F(a, 1 — a; 1; r^) 



4.10. Theorem. [22] Let < a,b < 1, max{a + 6, 1} < c < 1 + min{a,6}, and let 
Q be the quadrilateral in the upper half plane M. = {z E C : Im z > 0} with vertices 
0, 1, A and B, the interior angles at which are, respectively, bn, (c — b)TT, (1 — a)7r and 
(1 + a — c)7r. Then the conformal modulus of Q is given by 

(4.11) QM(A B, 0, 1) = M{Q) = X(r')/X(r), 
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where r G (0, 1) satisfies the equation 

^j,/2(c-a-6)^j^^ — a, c — 6;c+l — a — 6; r'^) 



(4-12) A-1 , , 



say, and 



^ E(c-&,1 -g) (fc^i,^).^ 
B{b,c-b) 



For a fixed complex number b with Im (6) > define the following function 
g{x,y) = QM(x + i ■ y,b, 0,1) for x G M, y > 0. This is well-defined only if the 
polygonal domain with vertices x + i ■ y, b, 0, 1 is positively oriented. This holds 
e.g. if Re (6) < and a: > 0. It is a natural question to study the level sets of 
the function g . This function tells us how the modulus of a polygonal quadrilateral 
changes when three vertices are kept fixed and the fourth one is moving. For instance, 
it was shown in [19] that the function decreases when we move the fourth vertex into 
certain directions. 

4.13. Trapezoid (Burnside [12]). In P pp. 237-239] so called square frame, the 
domain between two concentric squares with parallel sides, was considered. Such a 
domain can be split into 8 similar quadrilaterals, and we shall study here one such 
quadrilateral with vertices 1 + hi, {h — l)i, 0, and 1, h > 1. When h > 1 we have by 
[ini pp. 103-104], [H] 

(4.14) M(Q; 1 + hi, {h - l)i, 0, 1) = M{h) = X{r)/X{r') 

where 

ti — t2\'^ _i / TT \ _i /^C\ 

t^TTj ^ ^i = /^V2(^j' ^2 = /ir/4yj, c = 2h-l. 

Therefore, the quadrilateral can be conformally mapped onto the rectangle l+iM{h), 
iM{h), 0, 1, with the vertices corresponding to each other. It is clear that h — 1 < 



M{h) < h . The formula (4.12) has the following approximative version 

M{h) = h + c + 0(e~^^), c = -1/2 - log 2/tx ^ -0.720636 , 

given in [33]. As far as we know there is neither an explicit nor asymptotic formula 
for the case when the angle 7r/4 of the trapezoid is replaced by an angle equal to 
a G (0,7r/2). 

4.15. Numerical computation of elliptic integrals. The computation of the 
elliptic integrals is efficiently carried out by classical methods available in most pro- 
gramming environments. Numerical estimates for X(r), and hence for /Xi/2(r), are 
obtained very efficiently by the following recursive method. For r G (0, 1) let 

ao = l, bo = r'E (0,1), 

O'n+l = {0"n + ^n)/2, &n+l = V O'nbni 

Then the sequences (a„) and (6„) have the common limit 7r/(2lX(r)), and, for each 
y G (0, oo) we can approximate jjC^j^iv) numerically by the Newton- Raphson iteration. 
For details, see e.g. [4;, 3.22, 5.32] and [2?, 2.11]. 



ON MODULI OF RINGS AND QUADRILATERALS 



13 



5. Validation of algorithms: convex quadrilaterals 

Validation of the algorithms for the modulus of a quadrilateral will be discussed 
in two main cases: convex quadrilaterals and the case of a general polygonal quadri- 
lateral. In this section the case of a convex quadrilateral will be discussed for the 
following three algorithms: (a) the SC Toolbox in MATLAB written by DriscoU [T7] . 
(b) the AFEM software due to Samuelsson [9], (c) the hp-method of the present 
paper implemented in the Mathematica language using the double precision. The 
reference computation is carried out by the algorithm of |22], implemented in [22] 
in the Mathematica language (the algorithm QM[A,B] implementing the formula in 



Theorem 4.10). This implementation makes use of multiple precision arithmetic for 
root finding of a transcendental equation involving the hypergeometric function. All 
the SC Toolbox tests in this paper were carried out with the settings precision = 
le-14. 

5.1. Setup of the validation test. All our tests were carried out in the same 



fashion using the reciprocal identity (4.3) and considering a quadrilateral with the 



vertices a, b, 0, 1 with Ima > 0, Im6 > 0, and the line segments joining the vertices 
as the boundary arcs. The vertices b, 0, 1 were kept fixed and the vertex a varied over 
a rectangular region in the complex plane. The numerical value b = —0.2 + i- 1.2 was 
used and the lower left (upper right) corner of the rectangular region was 0.5 + i • 0.2 
(1.5 + 2-1.2). Examples of such quadrilaterals, along with some minimal meshes used 
in the computation, are illustrated in Figure |3} The test functional, based on the 



reciprocal identity (4.3), is 

(5.2) test(a, b) = |QM(a, b, 0, 1)QM((6 - l)/(a - 1), 1, 1/(1 - a), 0, l) - l| 

which vanishes identically. The values of this test functional are reported in Table 
[T] for the fixed value b = — 0.2 + Z-1.2 when a runs through the aforementioned 
rectangular region. A table of values of QM(m + in, i, 0,1), m,n = 1, . . . , 5 is given 
in [221 Tablel]. 



Table 1. Tests related to the convex quadrilateral, with (0.15,18)- 
meshes used in the hp-method. 



Method 


Error range 


AFEM 
SC Toolbox 
hp-method {p = 12) 
hp-method {p = 15) 
hp-method {p = 18) 


1.15 ■ 10-" 
1.11 ■ 10-^^ 
6.51 • 10-1^ 
2.22 • 10-^^ 
1.11 ■ 10-16 


1.41 ■ 10-« 
1.55 ■ 10-1^ 
7.84 ■ 10-9 

1.42 ■ 10-1° 
3.90-10-12 



5.3. The reference computation. We used the Mathematica script of 



for 



solving the equation in Theorem 4.10 for the computation of QM(a, b, 0, 1) in order 
to carry out the test. The conclusion was that the amplitude of the error was roughly 
10"!'' i.e. there was practically no error. Note that the quadrilateral here is not always 
convex. On the basis of numerical experiments, it seems that the reference method 
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Figure 6. Logarithm (with base 10) of errors over the domain 
[0.1, 2] X [0.1, 2], corresponding to values of p = 12, 15, 18 starting from 



above. The error estimate is obtained by using the identity (4.3). 



of |22] does also work in non-convex cases, but this has not been rigorously proved. 
These results were then compared to numerical values obtained by using AFEM, SC 
Toolbox and the hp-method. The results are summarized in Table [T] 

6. Validation: polygonal quadrilaterals 

In this section we will consider the validation of the algorithms for the modulus 
of a quadrilateral in the case of polygonal domains with g > 4 vertices. In the 
case considered in the previous section there was a reference computational method, 
providing the reference value for the moduli. There is no similar formula available 
for the general polygonal case. 

6.1. Setup of the validation test. All our tests were carried out in the same 



fashion as in the previous section, using the reciprocal identity (4.3). We selected a 
quadruple of points {zi, Z2, z^, z^} , which is a subset of the set of vertices defining 
the polygon D , and assume that these are positively oriented. Thus (D; zi, Z2, z^, z^j 



is a quadrilateral to which the reciprocal identity (4.3) applies. 



6.2. The notation cmodu(w, fci, /C2) and modu(w, fci, ^2) . Suppose that w is a 
vector of p complex numbers such that the points Wi, . . . , Wg, q > 5, are the vertices 
of a polygon D and that they define a positive orientation of the boundary. Choose 
indices ki,k2 G {1, . . . ,p — 1} with ki < k2 and set zi = Wk-^, Z2 = Wfe^+i, Z3 = Wk^, 
Z4 = Wk2+i ■ Then we define 

cmodu(w, ki, k2) = M(D; Zi, Z2, Z3, Z4) , modu(w, ki, ^2) = M{D; Z2, z^, z^, Zi) . 
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By the reciprocal relation (4.1) we have 

(6.3) cmodu(w, ki, ^2) ■ modu(w, /ci, = 1 . 

6.4. L-shaped region. The L-shaped region: 

L(a, 6, c,d) = LiU L2, Li = {z e C : < Re (z) < a, < Im (z) < b}, 

L2 = {z e C : <Re{z) < d, <lm{z) < c} , < d < a, < 6 < c , 

is a standard domain considered by several authors for various computational tasks. 
In the context of computation of the moduli it was investigated by Gaier [20j| and we 
will compare our results to his results. In the test cases all the vertices had integer 
coordinates in the range [1, 4] . Since we consider an integer coordinate domain, simple 
quadrilateral grid has the desired properties of the minimal mesh. An example of 
such a mesh is shown in Figure [9j The results are summarized in Table [2| and the 
potential functions are illustrated by Figure [7j 





Figure 7. Potential functions in the case of L-shaped region |6.4 
The vertices of the region Q are zi = (0,0), Z2 = (3,0), Z3 = 
(3,1), 24 = (2,1), Z5 = (2,2) and Zq = (0,2). Potential func- 
tions related to M{Q] Z2, z^, zq, z^) ^ 1.5081540958548603 (left), and 
M{Q; zi, Z2, Zi, ze) ^ 0.6630622181450123 (right), are illustrated. 



Figure 8. Minimal (see 3.16) meshes for domains of 6.4 and 7.1 
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Table 2. Tests of (6.3) for L-shaped regions (see 6.4 and Figure [?]), 
with (0.15, 12)-meshes used in the hp-method. 



Method 


Error range 


AFEM 
SC Toolbox 
hp-method {p = 12) 
hp-method {p = 16) 
hp-method {p = 20) 


1.80 ■ 10-1" 
2.22 • 10-1^ 
4.01 ■ 10-11 
8.03 ■ 10-13 
5.98 ■ 10-13 


7.10 ■ lO-i'^ 

2.58 ■ 10-1^ 

1.59 • 10-1° 
2.28 • 10-12 
1.80 ■ 10-12 



Table 3. Table for square in square, p = 16. 



a 


Capacity 


Error (hp) 


Error (SC) 


0.1 


2.83977741905223 


2.35 


10- 


'15 


3.17 


10- 


-14 


0.2 


4.134487024234081 


1.93 


10" 


-15 


2.10 


10" 


-15 


0.3 


5.632828000941654 


1.58 


10" 


-16 


2.69 


10" 


-16 


0.4 


7.5615315398105745 


1.17 


10- 


-15 


1.50 


10" 


-16 


0.5 


10.23409256936805 


1.74 


10- 


-16 


3.42 


10- 


-16 


0.6 


14.234879675824363 


7.49 


10- 


-16 


1.35 


10- 


-15 


0.7 


20.901581676413954 


0. 






2.89 


10" 


-16 


0.8 


34.23491519877346 


6.23 


10- 


-16 


3.61 


10" 


-16 


0.9 


74.23491519877882 


3.83 


10- 


-16 


8.31 


10- 


-15 



Table 4. Table for cross in square, p = 16. 



a 


b 


c 


Capacity (SC) 


Capacity (hp) 


Difference 


0.5 


1.2 


1.5 


21.94721953515564 


21.94721953515577 


5.99 


10" 


-15 


0.5 


1.0 


1.5 


14.00279904484107 


14.00279904484109 


8.88 


10" 


-16 


0.2 


0.7 


1.2 


9.186926595881523 


9.186926595881525 


1.93 


10- 


-16 


0.1 


0.8 


1.1 


11.256582318490887 


11.256582318490889 


1.58 


10- 


-16 


0.5 


0.6 


1.5 


7.323269585560689 


7.323269585567927 


9.88 


10- 


-13 


0.1 


1.2 


1.3 


23.13861453810508 


23.13861453810529 


8.91 


10- 


-15 



7. Ring domains 

In this section, we compare hp-FEM with exact values and with AFEM and SC 
Toolbox in certain ring domains. The square in square and cross in square cases were 
previously considered in ^9j and numerical values were reported in [9^ Tablel, Table 
4]. Our numerical results in Tables [3] and |4] provide 12 decimal places whereas in [S] 
only 6 decimal places were given. 

7.1. Square in square. We compute here the capacity of the ring domain with 
plates E = [-a, a] x [-a, a] and F = Coo \ ((-1,1) x (-1,1)), < a < 1. The 
results with SC and the hp-method with (0.15, 16)-meshes are summarized in Table 
[T] For computation of the capacity, the ring domain is first split into four similar 



quadrilaterals. For the potential function, see Figure 10 Note that in this case, the 



exact values of the potential are known, see (4.14) and the related trapezoid type 
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Figure 9. Meshing setup for cross in square. 

Figure 10. Potential functions: square in square and cross in square. 
Because of the symmetry, only one fourth of the picture is shown. 

quadrilateral example. Explicitly, with c = (1 — a)/(l + a) and 



u 



A'' 1/2 



the capacity equals 47r//ii/2( 



vr c 

T 

r) 



U — V 



U + V 



7.2. Cross in square. Let Gab = {{x,y) '■ \x\ < a, \y\ < b} U {{x,y) : |x| < b, \y\ < 
a}, and Gc = {{x,y) : |x| < c, \y\ < c}, where a < c and b < c. We compute the 
capacity of the ring domain R = Gc \ Gab- The results with SC and the hp-method 
with (0.15, 16)-meshes are summarized in Table |4j For computation of the capacity, 
the ring domain is again first split into four similar quadrilaterals. The mesh for the 



quadrilaterals is given in Figure [9| and the potential function is given in Figure [10 
The exact values are not known in this case. 

Since the underlying mesh topology remains constant in both examples above 
we have computed the results using exactly the same mesh template for every sub- 
problem, e.g. Figure |9] for Cross in square, a = 0.5,6 = 1.2, c = 1.5, except for the 
extremal cases in terms of element distortion a = 0.9 for the square in square, and 
the case a = 0.5, b = 0.6, c = 1.5 for cross in square. Thus, the results also measure 
the robustness of the method with respect to moderate element distortion. Also, in 
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both cases due to symmetry we have graded the mesh only to the reentrant corners 
of the domain. 

7.3. Rectangle in rectangle. Let Gabcd = {{x,y) '■ a ^ x < c,b < y < d} and 

G = {{x,y) :l<x<7, l<y< 4}. We compute the capacity of the ring domain 
R = G \ Gabcd- Here we consider a subset of possible cases when a, b,c,d G N. 
The results computed using the hp-method with (0.15, 16)-meshes are summarized 
in Tab le [s) The potential function for the case {a,b,c,d} = {4,1,6,2} is given in 
Figure 111 The exact values are not known in this case. 



Again, we have employed the same mesh template (simple quadrilateral grid as 
in Figure |9]) over the entire test set. Grading has been used in the corners of Gated 
only. From results of Table [5] we can also see that some of the configurations are 
symmetric in terms of capacity. In these cases the differences in the computed values 
is less than 10" 



-13 



Table 5. Table for rectangle in rectangle. 



a 


b 


c 


d 


P 


Capacity 


1 


1 


2 


2 


20 


5.210320385649294 


1 


1 


3 


2 


19 


6.746053277945276 


1 


1 


4 


2 


20 


8.27007839293125 


1 


1 


5 


2 


19 


9.86240917550835 


1 


1 


6 


2 


17 


11.89718127369752 


2 


1 


3 


2 


18 


4.692072335693745 


2 


1 


4 


2 


18 


6.232078709256309 


2 


1 


5 


2 


20 


7.827105378062926 


2 


1 


6 


2 


17 


9.86240917550835 


3 


1 


4 


2 


17 


4.621123827863167 


3 


1 


5 


2 


20 


6.232078709256313 


3 


1 


6 


2 


18 


8.2700783929313 


4 


1 


5 


2 


19 


4.69207233569376 


4 


1 


6 


2 


20 


6.746053277945233 


5 


1 


6 


2 


20 


5.210320385649318 
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8. Domains with curved boundaries 

In this section, we give further examples featuring domains with curved bound- 
aries. Simple examples of such domains are domains, where four or more points 
are connected with circular arcs. Some examples related to numerical methods and 
Schwarz-Christoffel formula for such domains can be found in the literature, e.g. [TT] . 
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4 



Figure 15. The quadrilateral (Qb; vr/12, 177r/12, 37r/2, 1), the mini- 
mal mesh and the potential function. The relative error is 5.17 ■ 10^^^. 




10-'" 



5.0 7.0 10.0 15.0 20.0 



5.0 7.0 10.0 15.0 20.0 



Figure 16. The p-convergence of the reciprocal error 
for the quadrilaterals (Qyi; 7r/12, 177r/12, 37r/2, 1) (left) and 
(Qb; vr/12, 177r/12, 37r/2, 1) (right). Logarithmic scale. 




Figure 17. Wave: the minimal mesh and the potential function for 8.10 



Our method has the advantage that even more general quadrilaterals can be consid- 
ered, as illustrated by examples given below. Here the meshing has been tuned by 
monitoring the rate of convergence in the polynomial degree. Both the minimal mesh 
and the scaling factor have been adjusted until exponential convergence in p has been 



observed. Two examples of exponential convergence, as predicted by Theorem 3.15 



are shown in Figure 16 
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Figure 18. Flower I: the minimal mesh and the potential function. 
See also 18.111 and Table El 





Figure 19. Flower II: the minimal mesh and the potential function. 
See also 18.111 and Table IH 

8.1. Circular quadrilaterals. The absolute ratio of four points a,b,c,d G C is 
defined as 

\a — c\\b — d] 



8.2 



|a, 6, c, d\ 



\a — 6| |c — d\ 

The main property of the absolute ratio is the Mobius invariance: 
(8.3) \a,b,c,d\ = \w{a),w{b),w{c),w{d)\, 

if w is a Mobius transformation 

kz + I 



14) 



w{z) 



{kn — ml ^ 0). 



mz + n 

Given zi, Z2, z^ on a circle (or on a line) and wi,W2, W3 on a circle (or on a line), there 



exists a Mobius transformation w with w(zi 



Wj, j = 1,2,3. 



8.5. Type A. Let us first consider a quadrilateral whose sides are circular arcs of 
intersecting orthogonal circles, i.e., angles are 7i/2. Let 0<a<b<c<2n and 
choose the points {1, e"^, e^, e'^} on the unit circle with the absolute ratio 

sin(6/2) sin((c-a)/2) 



.6) 



l,e^e^e^ 



u. 



sin(a/2) sin((c- 6)/2) 

Let Qa stand for the domain which is obtained from the unit disk by cutting away 
regions bounded by the two orthogonal arcs with endpoints {l,e*"} and {6*^,6*^^} , 
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respectively. Then Qa determines a quadrilateral [Qa] e°, e*, e'^, 1) . Using a suitable 



Mobius transformation and the invariance (8.3) we can map Qa onto the upper half 



of the annulus {2; G C : 1 < < t} and we obtain the following formula: 
(8.7) M(QA;e^e^e^l) = 7r/logt, 

i.e. a half of the modulus of the full annulus, where 

t = 2u-l + 2Vn2 - M, t > 1. 
The results are summarized in Table HI 



Table 6. Moduli of quadrilaterals [Qa] e 



imn/24: „i?i7r/24 pirn/24: 



1) for 



several integer triples (m, n, r) computed with the hp-method, p = 20. 



Nodes 


Reference 


Computed value 


Relat. error 


Recipr. error 


(2,10,12) 
(2,10,14) 
(4,12,18) 
(6,16,24) 
(8,22,32) 


0.7071508111121534 

0.8074514311467651 

1.0383251171675787 

1.170060906774661 

1.313262425617007 


0.7071508111121347 
0.8074514311467831 
1.0383251171675796 
1.1700609067746603 
1.3132624256170076 


2.64 ■ 10-14 
2.23 ■ 10-14 
8.55 ■ 10-16 
5.69 ■ 10-16 
5.07-10-16 


1.02 ■ 10-13 
2.55 ■ 10-14 
1.44 ■ 10-15 
2.22 ■ 10-15 
3.33 ■ 10-16 



8.8. Type B. Next we let the sides of the quadrilateral be circular arcs be of the 
unit disk, and in this case all the angles are equal to vr. Now the unit disk, together 
with the boundary points e"', e^, e^, 1 determines a quadrilateral which we denote by 
Q B ■ Using an auxiliary Mobius transformation of the unit disk onto the upper half 
plane we can readily express the modulus using the capacity of the Teichmiiller ring 
domain |4l Section 7] and express it as follows 

(8.9) M(gB;e^e^e^l) = It{u-1), 



where u is as in (8.6), and 

T{t) 



TV / fli /2{1 / VT+t) , 



t > 0, 



and fii/2ir) is as in (4.9), gives the conformal capacity of the plane Teichmiiller ring. 
The results are summarized in Table [7l 

Table 7. Moduh of quadrilaterals (Q^; e*™^/24^ g*n^/24^ 

several integer triples (m, n, r) computed with the hp-method, p = 20. 



Nodes 


Reference 


Computed value 


Relat. error 


Recipr. error 


(2,10,12) 
(2,10, 14) 
(4,12,18) 
(6,16,24) 
(8,22,32) 


0.5389714947317054 
0.5953434982171909 
0.7121629047455362 
0.7718690862645192 
0.8319009599091923 


0.5389714947624924 
0.5953434982359955 
0.7121629047457778 
0.7718690862646902 
0.8319009599093506 


5.71 ■ 10-11 
3.16 ■ 10-11 
3.39- 10-13 
2.22 ■ 10-13 
1.90 ■ 10-13 


7.68 ■ 10-11 
4.26 ■ 10-11 
6.06 • 10-13 
4.09 ■ 10-13 
3.48 ■ 10-13 



Next we consider non-convex examples featuring quadrilaterals with curved 
boundaries which are not circular segments. 
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8.10. Wave. Let Q = {{x,y) : < x < 1, sin(27rx)/4 < y < 1 + sin(27rx)/4}, and 
Zi = (0, 0), Z2 = (1, 0), ^3 = (1, 1), ^4 = (0, 1). Then the hp-method with p = 20 gives 
M{Q; Z2, z^, Z4, Zi) ~ 1.285385932609546. An error estimate based on the reciprocal 
identity (O) is 2.66 • 10"^^ 



For visuahzation, see Figure 17 



8.11. Flowers. Let Q be the domain bounded by the curve r{9) = 0.8 + tcos{mi9), 
< 6 < 27r, = 4, 6, 8 and t = 0.1 or t = 0.2 . Domains of this type are illustrated in 
Figures [Tsj and 19 We compute moduh of quadrilaterals M{Q; zi, Z2, Z3, z^j, where 
Zj = r{6j). We consider flower shaped quadrilaterals of type I with 6j = (j — l)7r/2 
for j = 1, 2, 3, 4, and type II, where 9i, 62, 9^ are as before, and 64 = 57r/4 (see Figures 
18 and 19). The numerical results are summarized in Tables Isl and |9j 



Table 8. Moduli of flower-shaped quadrilaterals {t, n) of type I com- 
puted with the /ip-method, p = 20. Note that, because of symmetry, 
it follows from (4.1) that the exact value of modulus is 1. 



n 


Error (t = 0.1) 


Error (t = 0.2) 


4 
6 
8 


3.18 ■ 10-14 
3.74- 10-" 
1.34- 10-13 


2.25 • 10-14 
8.45 ■ 10-11 
6.27- 10-11 



Table 9. Moduli of flower-shaped quadrilaterals (t, n) of type II com- 
puted with the hp-method, p = 20. The error estimate is obtained by 



using the reciprocal identity (4.1). 



n 


t 


Error 


Modulus 


4 


0.1 


2.00 


10" 


-15 


0.8196442147286799 


4 


0.2 


1.40 


10" 


-13 


0.8196441884805612 


6 


0.1 


2.34 


10" 


-14 


0.7896695654987764 


6 


0.2 


1.43 


10- 


-10 


0.7690460663235661 


8 


0.1 


9.05 


10" 


-14 


0.8196441884804566 


8 


0.2 


1.38 


10" 


-10 


0.8196441885295815 



9. Summary 

The computation of the moduli of quadrilaterals and ring domains with piece- 
wise smooth boundaries is a problem frequently occurring in various applications, 
see ^4j- There is no general method for such computations except perhaps the case 
of polygonal quadrilaterals when the SC Toolbox |T7l |18] may be considered as the 
"state-of-the-art" tool. For the case of ring domains there is no such general tool, but 
the adaptive finite element software AFEM of K. Samuelsson pj] has turned out to 
be effective in a number of cases reported in [9j. For the purposes of this paper the so 
called hp-FEM method implemented by H. Hakula, and first reported in this paper, 
is used in several examples with curvilinear boundaries where the previous methods 
do not apply. The hp-FEM method, applied to the computation of moduli of two 
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ring domains previously considered in [9] and reported in Tables |3] and |4| provide a 
significant improvement over the values reported in [9|. 

For experimental error estimate we have used so called reciprocal identity, which 
we have not seen used anywhere for the purpose of error estimation. It is our belief 
that this simple identity should be more widely known. It provides a criterion for 
estimating the error of numerical computation of the modulus of a quadrilateral for a 
very large class of simply connected domains, including those with curved boundaries. 
It seems that such a large class of examples has previously not been known for 
instance in the FEM community. These examples also enable one to experimentally 
demonstrate the theoretical convergence rates in nontrivial model problem cases as 
we have shown for the case of hp-FEM. 

For the very special case of convex polygons with four sides, the modulus of the 
corresponding quadrilateral is known as an analytic-numeric formula (4.11) by [22] 
and this is our starting point. We compare the performance of SC Toolbox, AFEM, 
and hp-FEM against the formula [22] and the reciprocal identity test. Next, again 
using SC Toolbox, AFEM, and hp-FEM, we consider polygonal quadrilaterals with 
more sides, L-shaped quadrilaterals and carry out similar comparision, using again 
the reciprocal identity as the test quantity. Thereafter, we discuss, now using AFEM, 
and hp-FEM, two classical cases of ring domains, the square frame and the cross in 
square ring domains previously considered e.g. in [9] where further references may be 
found. The error estimate in the square frame case is based on the well-known formula 
(4.14) whereas for the cross in square case we use SC Toolbox and the results from 
[H] as the comparision data. Finally, we also consider several cases of quadrilaterals 
with curvilinear boundaries, now only using the hp-FEM method, because the other 
methods mentioned above do not apply. 
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